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SYNOPSIS 


Dynamic stability of pipes ccnveying fluid has 
been studied using finite element analysis# The element 
matrices are obtained using GalerldLn’s method# It is 
shown that this method gives a systematic and straight 
forward way to account for natural boundary ccndi ticns 
unlike other methods. The results obtained through 
finite element analysis agree very well with those 
obtained by using classical methods. The stability of 
various pipe configurations, for example pipes with 
various types of elastic supports at both ends and 
fixed-pinned configuration, have been studied by this 
method# 


A very. general computer programme has been 
developed to solve these problems. An important 
feature of the programme is its flexibility. Solution 
for different boundary conditions can be obtained simply 
by changing the input cards. Also, generalization for 
ncn-uniform cross sections and/or multispan pipes with 
non-periodic and ncn-identical supports can be easily 
incorporated# 



CHAPTER 1 


OTRQDIICTION 

1*1.. Dynamic St a bility of Pipes Conveying 
i JFluid by gjM 

Lateral vibration of pipes conveying fluid has 
been studied extensively in the last two decades. The 
dynamic. behavior of these pipes subjected to transverse 
loading is of great importance in designing systems such 
as oil pipe lines, heat exchanger tubes, feed lines to 
-rocket motors and water turbines, etc. These problems 
belong to a branch of engineering mechanics which has come 
to be known as ’’flow-induced vibrations”. The first book 
on this subject has been published recently, Blevins [31* 

fluid flowing through a system is a source of 
energy that can induce vibrations and lead to instabilities. 
As the flow velocity increases, a certain value is reached 
at which the pipe looses stability by buckling and/or 
flutter, and the velocity at which this occurs first is 
called critical velocity. 

The stability of fluid conveying pipes is of 
practical importance because the natural frequency of a 
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pipe generally decreases with increasing velocity of fluid 
flow. The pipe may become susceptible to resonance if its 
natural frequency falls below certain limits. 

The analysis of such problems when done by classi- 
cal methods is quite involved even for uniform pipes with 
simple boundary conditions. Also for every case of the 
boundary condition a fresh characteristic equation has to 
be solved. 

For complex pipe line problems with non-uniform 
oross sections or multispan pipes with non-periodic, ncn- 
identical supports, finite- element analysis is found to be 
very convenient. Moreover, all types of boundary condi- 
tions can also be accommodated very easily. 

1 4 2 Review of Previous Work : 

A comprehensive review of the literature with 
extensive references is given in Paidoussis and Issid II 7 ] 
and Chen 1 5 1 * 

Historically, in the early 1 950’ s, observation 
of bending vibrations of the Trans - Arabian pipe-line 
prompted an investigation by Ashley and Haniland [11 . 

They studied the vibration of a simply supported pipe and 
suspected seme agen cy which damps the pipe line motion 
and tends to destroy its simple harmonic nature* The same 
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problem was studied independently by Housner {10] using a 
different approach. He found that at low fluid velocities 
the effect upon the vibrations of the pipe was negligible 
but at a certain high velocity the pipe buckled, like a 
column subjected to axial loading, A general study was 
made by Niordson [15 ] who used classical shell theory to 
deduce the equation of motion for axial, peripheral, and 
radial deformation of the shell, which led to the same 
conclusions regarding the stability of a pipe with simply 
supported ends. 

Long f 1 3 3 took the case of cantilever . pipes 
conveying fluid. He adapted an iterative procedure using 
a power series for the mode shape, which ^vas applicable 
to relatively small flow velocities and ecnfiimed experi- 
mentally that the forced motion of cantilever' l pipes were 
damped by internal flow in the range of flow velocities 
considered. 

Benjamin t 2 3 dealing with the general dynamical 
problem of articulated pipe system conveying fluid was the 
first to anticipate the phenomenon of unstable oscillations 
when such system possesses one free end. He produced a 
complete theory, supported by experiments for articulated 
pipe system. He showed that when the pipe is vertical 
both oscillation and buckling instabilities are possible 



in general, whereas when the motion is confined to a hori- 
zontal plane, that is when gravity is insignificant buckl- 
ing can not occur. Later Gregory and Paidoussis C 8] 
showed theoretically and experimentally that at sufficient- 
ly high flow velocities cantilever pipes are subjected to 
oscillatory instability only. Next Paidoussis 1 1 7 3 included 
the gravity forces and found that buckling instability is 
not possible in cases of hanging cantilevers only a flutter 
type instability is possible. This paradox was clarified 
by Paidoussis and Deksins [18] . 

Naguleswaran and Williams [16 1 studied the effect 
of internal pressure of the fluid theoretically and experi-, 
mentally. In the case of pipes supported at both ends 
the effect of internal pressure was found to be similar 
to that of flow velocity; that is, the pipe can buckle 
even at very low flow velocity by the action of sufficient- 
ly high internal pressure. 

More recently in 197*1 Chen [*+] studied the 
stability of a pipe conveying fluid, with the upstream 
end of the pipe fixed and the down stream end supported by 
a displacement spring, so that the boundary conditions 
were intermediate between clamped - free and clamped - 
pinned; accordingly both buckling and flutter were possible 
depending on the spring constant. 
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In a recent paper, Paidoussis and Issid [ d 9 ] 
considered the dynamics and stability of flexible pipes 
containing flowing fluid, where the flow velocity was 
either entirely constant, or had a small harmonic ccmpo- 
nent superposed. It was shorn that conservative systems 
were subjected not only to buckling but also to oscilla- 
tory instability (flutter) at higher flow velocities. The 
dynamics of a fluid- conveying pipe, clamped at the up- 
stream end and supported by a displacement and a rotational 
spring at the downstream end was studied by Lin and 
Ghen [ 12] , It was shown that the pipe may lose its 
stability by buckling, flutter, or both, depending cn the 
magnitudes of the displacement as well as rotational spring 
constants, 

1 * 3 Present Work, Objective and Scone 2 

The aim of the present work is to develop the 
finite element method for solving the problem of dynamics 
of pipes conveying fluid. 

In the second chapter the equations of motion (in 
the matrix form) are obtained by finite element analysis 
using Galerkin’s method. The stiffness, damping and mass 
matrices for a typical element are obtained. Here it is 
seen that Galerkin l s method. of finite element analysis 



6 


gives very distinctly a way to account for natural boundary 
conditions. It is noted here that the stiffness and the 
mass matrices are symmetric but the damping matrix is n cn- 
symmetric. These matrices are arranged to form the stan- 
dard dynamical matrix, Meirovitch [14 ] . 

Using a library subroutine the complex eigenvalues 
of the dynamical matrix are obtained for various boundary 
conditions of the pipe problem. These results are given 
in the third chapter and discussed. The mechanism of 
instability is also discussed at the end of this chapter. 

In the fourth chapter the conclusions are reported. 
Lastly the listing of the detailed Computer Programme is 
presented .in the appendix. 



CHAPTER 2 


THEORY AND METHOD OF SOLUTION 

In this chapter the fourth order partial diffe- 
rential equation governing the motion of the pipe con- 
veying fluid at velocity V is cast into simultaneous 
ordinary differential equations by finite element analy- 
sis using Galerkin’s method. These equations are further 
rearranged into the standard eigenvalue problem and the 
complex eigenvalues obtained. 

Finite ELement Analysis ; 

Consider a horizontal pipe conveying fluid at 
a velocity Y* The equation of motion for the small 
transverse vibration, can be found in Paidoussis and 
Issid [ I 9 ] . Neglecting internal damping and gravity it 
becomes 


(m p + + 2m f V + EE ^ 


+ (m f V 2 - T + P A f ) = 0 


( 2 . 1 ) 


8X 


and the boundary conditions are, 
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3 ' 

Tr«"r 3 Z m jbZ 

11 L? - T h 


+ k 


dl 


Oor z = 0 at x = 0 


El 


3 

3X 2 


- k t i 


dz 

3x 


= 0 or 


9 z 
3 x 


= 0 at x = 0 


3 

KE — - T r" - k ~ z = 0 or a = 0 at x = L 
3x 3 3 x d2 

2 

EI + kj. 2 = 0 or |^ = 0 at x = L 

t » « (2 * 2 ) 


where 

nif = Mass of the fluid per unit length of the pipe* 

m = Mass of the pipe per unit length* 

ir 

V = Velocity of the fluid passing through the pipe* 

A^ = Area of the passage through which the fluid passes. 

T = Externally applied tension to the pipe* 

P = Pressure of the fluid. 

E = Modulus of elasticity of the pipe material. 

I = Moment of inertia of the pipe. 

L = Length of the pipe, 
z = Transverse displacement of the pipe, 
x = Axial co-ordinate of the pipe. 


k 41 i 
*c(2 


Displacement spring constants at left and right 
ends respectively. 


k t1 » 
k t2 


= Rotational spring constants at left and right 
ends respectively. 



9 


It is desired to solve this problem by finite- 
element analysis using Galerkin's method* 


Let the beam have m nodes with m-1 elements 
as shown in Fig, 1(b). The typical finite element is 
shown in Fig.1 (c). 


The transverse displacement z within the ele- 
ment is a function of space and time. It is assumed that 
over the element z varies as 


(g\ f 

(x,t) = %(x) z,(t) = j N(x) i { z(t) } 

i=1 1 1 ,x 

(2.3) 

where f is the number of degrees of freedom assigned 
to the element (e) and z^(t) are the discrete values at 
the nodes. 


Applying Galerkin * s critericn^as done by Szabo 
and Lee [20 ] and Huebner [11 ] , we get. 



(e) 


+ 


3 i 2 z (e ) 

3 x 3t 



(e) 


] dx = 0 


i 1 ,2 ,f (2.4) 


The second, third and fourth terms of Eq.n. (2.4) 
are integrated by parts to introduce natural boundary 
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ccndticns. At the same time this leads to less stringent 
requirements for the interpolation function [_Nj . 

The second term gives, 


1 

0 


2 (e) 


1 


1 3 N ± 


0 


3 X 


~~ ) dx = 

3 X 3 t 

N . 2 Mf V if- 

! 0 

a -7 ( e ) 

2m f v tt 

dx « • * 

(2.5) 


The third term of equation (2 « 5 +) is integrated twice by 
parts, to get, 

1 (Kt 3lfz<e) ' ' 33z(e) 11 

0 


/ N, (El — jf ) dx = N. El ~ 

n 1 1 3 t3 


3 x 


3N-? o ■* n «2. 


0 


l-L 2 *™ I 1 +f 19 n i , T 3 2 ^ e > „ 

+ I — o — o dx 

0 0 3 x^ 3 xr 


ee 


( 2 . 6 ) 


The fourth term yields, 

d o a 2 (e) 

/ N- (m f IT - T + P ju) 2-3 dx 
0 x 3x 


= % (m f V 2 _ T + P Af ) 


lz (e) I 1 




0 


1 31, 2 (s) 

■ ; 0 3# (m f 7 - 1 +P %•) If- (2,7) 


Substituting Eqns. (2.5), (2.6) and (2.7) into Eqn. (2*4), 
one gets, 
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r 1 r/ 9 2 7 (^) ^-( 0 ) 

I {(m p + m f ) IT,. 1-| -a.j.VjJ i| 


‘P 


+ El 


3t 

a 2 % aVe) 

P p 

3z 3X^ 


f 3 x 3 t 


1 , o m x 3 N. (e) 

- (m^ TT- T + P Af) — - ~ } 
1 1 jz 3x 


dx 


, a 2 (e) 1 a3 z (e)il 

N, (2m V) |f | + N, El 

1 f 3t ’O 1 3x 3 i 0 

8H ± a 2 (e) i 1 p a (e),l 

— ^ El -M i + N ± (m f V- T + P Af)|| =0 


3X 


f o 

i = 1 ,2, ,f (2.8) 

Substituting the value of from Eqn. (2.3) in Eqn,(2.8), 

one gets, 


1 a 2„ ne- 

/ [ (m p + m f ) { N } [NJ { L* > 


3f 


n © 

-2m f V t »' } 1_HT J {||} 


+ El { If"} !_N" f { zi 


ne 


2 Ti© 

- (m f V - T + P Af) { N f } jjr'j (z) 1 dx 


- (N } El 


a3 z (e) f l 


3X- 


0 


HO (m f l / 2 - T + P Ap) 

T 3x i 

. (e),l , .2 ,1 

{ U } ( 2 m^. V) ^ | + { E } El — ^ | 


(e) ,1 
0 


3x '0 

( 2 . 9 ) 
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Eqn. ( 2 . 9 ) can be written in the following fonn, 

r , (e ) .. ne (e) , n.e (e) ne 

Cm] { z > + [ c ] { z } + [k] { z } 


El 


3„(e) ,1 


8 J z 




! IL EI ~ 

* -r-J 


3x3 ' 0 } 

3*(e) ,1 I 

!o 


\ 


| N f EI E_z 
! 3 x J 


9x 


ji.(e) ,1 r 


0 




" N -| (m f "T 2 - T + P A f )-| 


z 
3 x 


(e) .1 


i 

< 


N n (m, -V 2 - T + 


Af)- 


_3_z 

3 X 


(e) 


N f (m f V 2 - T + P a*.)JL2 


(e) 


f 3 x 




“1 

^ (e ) . 1 

N 2 (2m f V) 3 z » 


0 

1 

0 

1 

0 






/%-! rrr 

K* EE p 
3 x 


3 2 zC- e ^ I 1 "I 


3 t 


r 0 


? 


2 

<n' EE - — | 

1 a 3 X 2 ' 0 


*0 ! 

. 1 


*f (2 “f H 


(e), 1 


0 


. a 2 (e) , 1 

1 Hi- ei 3 _a' j 
L 3 X 2 1 0 J 


(2«1q) 
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where 

= Mass matrix of the element 

= / 1 (in +m r ) { N } In! dx (2.11) 

X=Q p x ' 

= Damping matrix of the element 

1 ‘ f 

= / - 2m f V { IT } j N ! dx (2.12) 

x=0 ' 

= Stiffness matrix of the element 

= J 1 [El {if} jjf]- (m f V 2 ~T-*P N' h dx 

(2.13) 

According to Szabo and Lee [20 3 the column 
matrices, on the right hand side of the Eqns- (2.1o), are 
the measures of the error of discretization of the finite- 
element method. But it is these matrices that introduce 
the natural boundary conditions of the problem and as such 
must be taken into consideration. This fact was recog- 
nized by Huebner till , but there is still sane confusion 
when elasticity problems are considered. 

While choosing the shape function [_N J the follow- 
ing conditions must be met to ensure mcnotcnic convergence, 
Huebner [11 1 , All uniform states of z and its partial 
derivatives upto the highest order appearing in Eqn. (2,8) 


[m] 


(e) 


[c] 


(e) 


[kj 


(e) 
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should have a representation in z 1 ^ when in the limit 
the element size shrinks to zero. Also at the elemental 
interfaces z and any of its partial derivative cne 
order less than the highest appearing in integrants of 
Eqn. (2.8) must be continuous. Elements satisfying the 
first requirement are known as complete and elements 
satisfying the second requirement are known as compatible 
in the literature. 


The following shape function commonly used for 
beam elements has been taken for this problem as it satis- 
fies the above requirements, Desai and Able [ 6 1 


H 

U? (3 - 2 50 

4 5 2* 

? 2 (3 - 2 

5 2) - 5 1 5 2<J 





a - * 4 

(2-14) 

where 


= 1 I 

\ 1 ~ 1 > 

5 = x 

2 1 



and 

1 

- length of the beam 

element 


Substituting the shape function (2.14) into equations (2.11), 

(2.12) 

and ( 2 . 13 ) and integrating, 

cne gets, 



(e) (m-p +m f ) 

1 

X 





1 56 22 1 

54 

- 1 3 1 i 

1 



4 l 2 

13 1 

156 

- 3 l 2 

- 22 1 

! 

E 

(2.15) 



Symmetric 


4 l 2 ; 

1 

L 



i im _ 


i 

I 



1 ? 


and 


i 

r 30 

6 1 

+ 30 

- 61 i 


1-61 

0 

6 l 

- I 2 J 


| - 30 - 

6 1 

- 30 

6 1 ! 


6 l 

l 2 

- 6 1 

o 1 

— 






* * * 

(2.16) 

i 12 

6 1 

- 12 

6' 1 ~i 


If l 2 - 6 1 2 l 2 


Ek] 


(e) _ ei 
l 3 


I - P Af - m f v 

+ ( 2 ± \ 

^ 301 j 


12 

Symmetric 
36 3 1 - 36 

4 l 2 - 3 1 
36 • 

Symmetric 


- 6 1 
4 l 2 

3 1 
-l 2 
-3 1 
4 i 2 


1 


(2.17) 


and Eqns. (2,1o) become 


(e) 


ne 


(e) ne 

* ' r • i 


(e) 


ne 


-"EE’ 


3„(e) 


=7 


u 


[m] v '{z} +[c] (z> +[k] (z> = -< 


-(m f V 2 - T +P Ap) 


+EE- 


© 


0 


< 


+(m f ? - T + P If ) 6 

0 


J 


k 


J 




-2m f V z. 

-L J 

0 

2mf> V z‘ k 
0 


-EE 2 




El — I 

s- 

(2.18) 


0 

V e > 

x 3 

' ° n 

i <e) ] 

X 2 Ij V 
o 

2 ~( e ) 


L 


k 
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For i}fte Whole domain of the problem, Eqns. (2.1 8) , becane 
|MJ BB {V> n + [ c ] BB {2} n +(KJ bb < Z ) n « 




where [M]g B ? [ S] BB ? t K ]g B are the assembled mass j 
damping and stiffness matrices, respectively, before the 


boundary conditions are applied. For computational ease, 
Eqns# (2.19) are rearranged as fellows; 
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It should he noted that 2m f V z-j and 2m f V z m wi II go 
to Qj 1 andC^^ respectively, whereas 


/ 
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(m f V 2 - T + P Af ) 0 1 and (m f V 2 - T + P A f ) Q m will go 
to K-j ^2 and K 2 m _i ? 2m respectively. 

After applying boundary condition Eqn. (2.2o) 
can be written in the form, 

[M ] {z > n + [ C ] { z } n + [K ] { z } n = {O'} (2.21) 

2.2 Method of Solution : 

The solution of these homogeneous set of diffe- 
rential equations is obtained as detailed in Meirovitch |~14-J . 
It should be noted that it is not necessary, as Meiro- 
vitch [ 14-] requires that [ M 1 , t C ] , t K 1 matrices be 
symmetric, see Frazer, Duncan and Collar [7 ] , In the 
following, Eqns. (2,21) are taken as n equations. 

Using generalised velocities ( z } as auxiliary 
variables, n second order ordinary differential equations 
(2.21) are converted to a set of 2n first order ordinary 
differential equations, Meirovitch [ 1 4-1 

[M] (7 (t) } + [K] { y(t)> = { Y(t) } (2.22) 

where 

f{ z (t)}] 

{ y(t) } ={ }snd { T(t) }= {0} ( 2 . 23 ) 

(*)}{ 

are column matrices consisting of 2n elements representing 
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generalized, coordinates and. generalized forces respectively, 
and 


[M ] = 


[0] [ M] 
CM] [G ] 


and 


[K ] 


are real matrices of order 2n. 



[ 0 ] 

[K] 


( 2 , 2 )+) 


Now we have the differential equation which in the 
matrix form is , 

CM] {y(t)> + C K] { y(t) } = (0} (2,25) 

Let {y(t ) } = e at {$ } (2.26) 

where fi represents a vector consisting of 2n cons- 
tant elements. Substitution of Sqn. (2.26) in Eqns.(2.25) 
leads to the eigenvalue problem. 


a CM H 

+ Ck 3 (16 > 

= (o> 



(2.27) 

which can be written as 





CD! {$} 

= ~ ijt 5 ) 

a 




(2.28) 

where 







”•[0] 

II] 


i 

CD] = - 

CK ] t M ] = 

-1 






- CK] t M] 

- CK] 

tc] 

\ 





9 

( 2 . 29 ) 

This plays 

the role of the dynamical matrix. 
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The eigenvalues of the above dynamical matrix 
(2.29) will in general be complex. Their real parts 
determine the stability of the system and the imaginary- 
parts give the frequency of vibration. If the real part 
is negative, the system will be stable. The system will 
have oscilatory instability (flutter) If real part Is 
positive and imaginary part is non-zero- On the other 
hand if this imaginary part is zero, the system will have 
buckling instability. 

In this formulation the effects of various types 
of geometrical and natural boundary conditions as well as 
the effect of axial load, fluid pressure and velocity cn 
the pipe stability can be studied very easily. 

1 computer programme for calculating the complex 
eigenvalue and eigenvectors of a dynamical matrix was 
written, Meirovitch [I 1 !-] . Later cn an in-built library 

subroutine was used to calculate the complex eigenvalues. 



CHAPTER 3 


RESULTS AND DISCUSSION 

This chapter deals with the numerical results 
obtained for pipes conveying fluid with different boun- 
dary conditions* Complex frequencies of the lowest four 
modes have been found as a function of the velocity of 
the fluid flow through the pipe. Fran the analytical 
point of view the mechanism of instability has been dis- 
cussed. The results are given in terms of the dimension- 
less parameters usually used in the literature. 

3* "l Dimensionless Parameters ; 

mf 

B = Dimensionless mass = ~ - + — 

mp + mp 

u = Dimensionless fluid flow velocity = 

y EI/m f 

Re(w)= Real part of dimensionless frequency 



+m p ) I>) 

Im(u)= Imaginary part of dimensionless frequency 

yj ■ El/ ( (mf +m p ) I>) 
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k d 

a d = Dimensionless displacement spring factor = — 

ICj. L 

a t = Dimensionless rotational spring factor = 

El 

.3* 2 Number of Finite Elements : 

To begin with one must decide the number of 
finite elements to be used in the analysis. First, for 
every boundary condition results were obtained for a few 
velocities with the number of elements varying between 
3 and 7 . It was found that after five elements there was 
hardly any change in the results in all the cases. The- shcfple 
results for pinned-pinned and cantilever pipe are shown in 
Table 1. All results reported in this dissertations are 
for the pipe with five finite elements. 

3l.3- GsclelIss Fre qu enc ie s , op ines : 

The calculated frequencies of the pipe problem 
for different boundary conditions are plotted in Argand 
diagrams. The imaginary part of the eigenvalue 3m(w ) 
which gives the frequency of the vibration is plotted cn 
the abcissa and the real part Re( ") is plotted cn the 
ordinate. The various parameters are u, $ , and a 

3 . 3-1 Pinned-Pinned Pipe 

Figure ( 2 ) shows the dimensionless canplex fre- 
quencies of a pinned-pinned pipe with dimensionless flow 




.814- 0 35*925 9*699 0 o 1 1 8.924 

.793 0 35*466 9.383 0 0 17.757 

.78 7 0 35.329 9*285 0 0 17*4-04- 


2.785 0 35*279 9*247 0 0 17.269 

2 .7 £4- 0 35.256 9.231 0 0 17.210 


23 



vO t>~ 



Number of Elements Vs Frequencies 
(b) Cantilever Pipe ( 6 = 0*5) 
















2 ? 


velocity u as a parameter for g = 0 . 5 . It is noted that 
the results match with those obtained by the classical 
method, Paidcussis and Issid [ I9 ] . Conplete agreement 
was observed for the case 3 = 0.2 also, though it is 
not presented here* 

3«3«2 Fixed-Fixed pipe 

Figure (3) shows the dimensionless complex fre- 
quencies of a fixed- fixed pipe with dimensionless flow- 
velocity u as a parameter for 3 = Ch 5* These results 
also'match with those by classical methods, Paidoussis 
and Issid t 19 1 . 

3 * 3 • 3 Can til ev er Pipe 

Figure (4-) shows the dimensionless complex fre- 
quencies for a cantilever pipe with dimensionless flow- 
velocity u as a parameter for 6 = 0.2. Here again our 
results agree with those obtained by the classical method, 
Gregory and Paidoussis [8] . The results for the canti- 
lever pipe for 3 = 0.295 also match with the classical 
one in the same reference. 


Here it is very Important to note that the cdumn 
matrices ' { Ni 2m t 7 |q > and { (m f V 2 - T + P Af)£2 

found in Eqn. (2.8) must be taken into account in damping 


(e) 


0 
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and stiffhess matrices, respectively. In the analysis 
these can not be neglected as suggested in the work of 
Szabo and Lee [20 ] . Results obtained neglecting these 
terms are shown in Fig (5)* This figure clearly indi- 
cates the large difference from the true results, 
(compare Fig. (!+)). 



The boundary conditions for this case are 

z = 0 and = 0 ? a-t x — 0 

0 (3.1) 

3 2 
El | _ ku?* = 0 and El £ = o ; at x = L 
3 X 3 3 'xr 

Figure (5) shows the dimensionless complex fre- 
quencies for the above mentioned system with dimensionless 
flow velocity u as a parameter for 0= 0.6 and a ^ = 1oO. 
From the figure it is clear that the results are the same 
as those of the classical method, Chen € ^ 3 . Results 

were also found to match for the case 6= 0.2 and = 10« 

Here also a large difference in the results was 

r 2 2^) I 1 

observed when the column matrices {Hi 2m f V — r } 

/ v 2. ^ u 

and {Ni(m f V 2 - T + P %) || S I Q }™ere not added in 
damping and stiffhess matrices. 
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3 « 3 • 5 E.i.Pe with cne End. Fixed and the other 
Supported by a Displacement and a 
Rotational Spring 

Boundary conditions for this case are 


z - 0 and <L2 = q ; at x = 0 

3 X 

3 - 2 

EE 3 — | - z = 0 and El 3 — § + k t2 ~| = 0 ; at x = L 

a x3 3^2 

• • • (3*2) 

Figure (7) shows the dimensionless complex fre- 
quencies for the system mentioned above with dimension- 
less flew velocity u as a parameter for 3 = 0.6, a ^ = 50 , 
a £ = 25 , and it is found that the results obtained by 
finite element method match with those by the classical 
method, Lin and Chen [12 ] . Agreement was also observed 
for the case of 6 = 0*2, = 1o and “ ^ = 5* 


{ Ni 2m f V 3 ~§ (e) 

O u 


It is to be noted again that column matrices 

T_ 

q 3 and ( Ni (mp V 2 _ 1 + p Jt^) °-£ 


3 „(e) 


3 x 


1 

0 


have to be added to damping and stiffness matrices, other- 
wise large discrepencies are found in the results* 


3*3*6 Pipe with One Bid Fixed and the Other Pinned 
Boundary conditions for this case are 

3 7 

z = 0 and -£ = 0 : at x = 0 

9 x 9 

9 2 z 

z = 0 and El = 0 
9x 2 


; at x = L 


(3*3) 
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Figure ( 8 ) shows the dimensicnless complex fre- 
quencies with dimensicnless flow velocity u as a para- 
meter for 3= 0,5* In the figure it is to be noted that, 
with increasing flow velocity, the frequency of the sys- 
tem decreases and the frequency of the first mode becomes 
zero at u =4-, 5 which is the first critical flow velocity 
for buckling. But at u *7-7 the system regains sta- 
bility in the first mode. At slightly higher flow velo- 
city at us 7*9 the first and the second mode loci 
coalesce cn the t Im( w) ] - axis, and coupled - mode 
flutter with first and second modes occurs. As the flow 
velocity u is further increased the frequency of vibra- 
tion vanishes at u =11,1 for the first mode. This is 
the second critical velocity (budding). At a slightly 
higher velocity u * 11. 2 the coupled mode flutter occurs 
in the second and third modes# 

In the case of this boundary condition it was 
noted that the critical flow velocity does not depend 
upon fj as in the cp.se of pinned-pinned pipes and fixed- 
fixed pipes, Paidous^is and Issidf 19 ] . Here the cri- 

tical velocity is 4.£ and is associated with buckling. 

It is interesting to note that this value lies between 
those of pinned-pinned ( tt ) and fixed-fixed pipes (2 ir ), 
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1*3*7 -JPiPe with Both the Ends Supported by a 
displacement and a Rotational Spring 
(Symmetric Boundaries*) 

Boundary conditions for this case are 



a 3 z 

+ k ai z = 0 


J2. 




El 

and 

El 

9 x 2 

- k ti 

az = 0 
ax 

; at x = 0 

El 

3 3* ' 

— % “ ^(32 z * 0 
9X° 

and 

3x 2 

+ k t 2 

a_z = 0 
ax 

; at x = L 


. • (3-4) 

Figure ( 9 ) shows, for B = 0 - 2 , = 10 and 

= 5 > the dimensionless complex frequencies for the 
above system with dimensionless flow velocity u as a 
parameter. With the increase in flow velocity u the 
frequency decreases and becomes zero at u « 4.25 in the 
first mode, which is the first critical flow velocity for 
buckling. At flow velocity u -4.8 the frequency be- 
comes zero in the second mode. As the flow velocity is 
increased further, the loci of the first and second modes 
coalesce on the [ Be( w ) ] - axis and leave the axis at 

symmetrical points, representing the coupled-mode flutter 
with the first and second modes. 

3. 3« 8 Pipe with One End Supported by a Displacement 
Soring and a Rotational Spring and the' Other 
Bid by a Displacement Spring Only (Ncn- 
Svmmetric Boundaries ) 


Boundary conditions for this case are 
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EI 3 x 3 + kd1 Z 
^3 - k a 2 z 


0 and EI- 


3 2 z 


ax 2 

0 and S^l 
9X 2 


Ir 3 Z 

K ti ax 0 


at x = 0 


0 ; at x - L 

( 3 . 5 ) 


Figure (10) shows, for 3 = 0.2, a d = 10 and 
a t ~ 5 j the dimensionless complex frequencies with dimen- 
sionless flow velocity u as a parameter. As usual in 
the beginning the frequency of vibration decreases with 


increase in fluid flow velocity. At u= 4, the pipe 
becomes unstable due to flutter in the second mode and 
is the first critical velocity for this case- With further 
increase in the flow velocity at u - b.2^ the frequency 
become zero in the first mode but still the system remains 
unstable in the second mode with flutter. The first and 
third modes always remain stable but the second mode is 
still unstable (flutter), moreover at flow velocity 
u = 10.2. there occurs flutter instability in the fourth 
mode as well. 


!• It 9 - EiPe with One, Sid Supported by a Displacement 

and a Rotational Spring and the Other EhdL by a 
Hotatlcnal Spring only C Mon -Symmetric Boundaries ') 

Boundary conditions for this case are 


EI § — | + z = 0 and 

EI £ _ jj. 

3X 2 « 

3-2 = o 

3 x 

; at x = 0 

* 3 2 



EE 3 — 5 =o and ElM 
^x- 3 3r 

+ k. 2 ^ =o 


; at x = L 



• « ' * 

(3.6) 
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Figure (11) shows for £=0.2, a ^~l0 and. 

a £ = 5 the dimensionless complex frequencies with dimen- 
sionless flow velocity u as a parameter. It is seen 
that this system is always unstable (flutter) in the 
second mode, 

3,4- The Mechanism of Instability : 

The mechanism of buckling and flutter instabili- 
ties of a pipe conveying fluid can be explained as 
follows: 


1*4.1 Buckling Instabilities 

It can be seen that the term (nif V 2 _ T + P Af)^-§ 

dx 

in the governing equation of motion act as an equivalent 
compressive force* As the fluid flow velocity increases, 
this compressive force increase thus reducing natural 
frequencies of the system. It may become zero and may 
lead to the buckling instability. Critical flow velocity 
may be found by the application of Euler’s method of 
equilibrium, Paidoussis and Issid [I 9 ] . But these 

results based on a static analysis may not hold true for 
the dynamic case ? Paidoussis and Issid t^]. 

Physically one can also visualize m-p V 2 § as 
the centrifugal force. When this force increases due to 
flow velocity and dominates over the flexural restoring 
force, the pipe buckles. 


CkNfb ; 

I to f y No. J" V 
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3« 1 +'*2 Flutter Instability 

This can be explained by energy considerations. 
When the energy transfer fran fluid to the pipe and vice 
versa exactly balances in one complete cycle of oscilla- 
tion, the condition is one of neutral instability, that 
is dynamic equilibrium. But when the former is more than 
the latter, then the amplitude of oscillation increases 
without limit* But when the energy transfer by the pipe 
to the fluid is more than the energy transfer from the 
fluid to the pipe, the oscillations are damped* 

.Analytically flutter may be attributed to the 
Coriolis component in the governing equation of the sys- 
tem. If the frequencies of two modes become equal, the 
flutter is called a coupled - mode flutter and oscilla- 
tions are of the form at exp (iwt), Paidoussis and 
Issid[ 19 ] . So far these forms of oscillations were 
observed in conservative systems alone. Now it is observed 
that they can exist in non- conservative systems as well. 

For example sec Fig. (9) for pipes with symmetrical spring 
supports. 



CHAPTER If 


CONCLUSIONS 


Based on the analysis of chapter 3 , I have 
reached the following conclusions. 


(1) 


( 2 ) 


The results of pipe problems by the finite 
element method are found to match with the 
results obtained by the classical method for 
the cases of pinned- pinned, fixed-fixed, canti- 
lever, cantilever with free end supported by 
displacement spring and cantilever with the 
free end supported by displacement and rotational 
springs. This shows that finite element equations 
of motion developed are correct - 


, 3 z ( e ) 

The column matrices t N n - V r 

1 


1 

0 


} and 


0 


} , found after 


% (m f V 2 - T +P A f )^ (e) 
applying Galerkin’s Method to the governing 
differential equation of motion, must be included 
in the damping and stiffness matrices respec- 
tively. These column matrices can not be neg- 
lected as is suggested by the work of Szabo and 



The stability of the pipe with its upstream and 
fixed and the downstream end pinned has been 
studied. It is found that the system looses 
stability by buckling at a velocity u A. 5 in 
the first mode and this critical velocity is 
independent of mass ratio B , as in the case of 
pinned- pinned pipes. Flutter is also possible 
at higher velocities in this case also. 

From the stability point of view this system 
should be designed such that the velocity of 
the fluid passing through the pipe does not 
exceed u = b . 5”. 

It is shown in Figures (9, 10, 11) that the pipe 
may loose its stability by buckling f flutter or 
both (coupled-mode), depending upon the magnitude 
of the displacement as well as the rotational 
spring constants.. One thing is interesting to 
note that for the case when the spring supports 
are symmetric the instability occures in buckling 
first, whereas in the case of non- symmetric spring 
supports the likelyhood of instability is by 
flutter in the second mode. 



35 


(5) The finite element method for solving pipe prob- 

lems has been found to be very accurate. By this 
method problems with different boundary conditions 
can be solved very easily. The results of the 
desired problem can be obtained just by changing 
the input data of the computer programme. 
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THE COMPUTER PROG HAM 

The computer programme gives the complex eigen- 
values of the assembled dynamical matrix. The programme 
’was written in Fortran I V and run on an 3BH 70’-^+ computer 
with 32K words of core storage. The programme is quite 
general except that same statements will be required to 
be changed for different boundary conditions. 

A Fortran listing of the programme is given in 
Appendix -II . 

T.l Description of Programme ; 

The MAIN programme reads and prints various 
data and outputs results. Other subroutines assemble 
the matrices and after applying boundary conditions it 
forms the dynamical matrix. The complex eigenvalues are 
calculated by calling a library subroutine 

Sub routine Name Description 

ASMEL Called by MAIN, assembles the ele- 

mental mass, damping and stiffness 
matrices in turn. 
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Description 

Called "by MAIN, applies the boundary 
conditions to the assembled mass, 
damping and stiffness matrices 
Called by MAIN, forms the dynamical 
matrix of the problem 
Called by DYNMAT, gives matrix inverse 
Called by MAIN, determines eigenvalues 
and eigenvectors of a general real 
matrix. The binary version of this 
programme is stored cn 1 301 disk file. 
The name of the file is ZINDLA. For 
further details refer to Communication 
of ACM Journal, Vol 2, I 968 , P. 82o 

LIST OF PRINCIPAL VARIABLES 

FORTRAN Implementation ; 

Program Symbol Definition 

(MAIN) 

ALPHA Dimensionless displacement spring 

constant 

BEL Length of the beam 

BETA Dimensionless mass 


Subroutine Name 
BCN DRY 

DYNliAT 

MAT IN V 
EIGENP 
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Program Symbol 


Definition 


C0MNC2 

CQMNIC3 


C0MKK4- 

C0MNK6 


CCMNK7 

C0MNK8 

DI 

DIM 

DOUT 

E 

ECM 

EKM 

H4M 

FACHON 

gama 


Factor which adds to damping matrix 
(see Eqn. 2.2o) 

Displacement spring factor at left 
end 

Rotational spring factor at left end 
Displacement spring factor at right 
end 

Rotational spring factor at right 
end 

Factor which adds to stiffness matrix 
(see Eqn. 2.2o) 

Moment of inertia of pipe 
Internal diameter of pipe 
External diameter of pipe 
Modulus of elasticity of the pipe 
material 

Elemental damping matrix C+,4) 
Elemental stiffness matrix 0+,*+) 
Elemental mass matrix 0+,*f) 

Dividing factor which gives dimension 
less eigenvalues 
Dimensionless rotational spring 
constant 
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Program Symbol 

Definition 

KB 

Factor gives the No. of constraints 

on the boundary 

NEL 

Number of finite elements of the beam 

NKB 

Matrix size after applying boundary 

condition 

m 

Total number of degrees of freedom 

p 

Fluid pressure 

HOF 

Density of the fluid 

ROT 

Density of the pipe material 

STIFAX 

Stiffness of the displacement spring 

STIFTO 

Stiffness of the rotational spring 

V 

Velocity of the fluid flow 

VN or 

Non-dimensional fluid flow velocity 

SUBROUTINE ASAMBL 

AM 

Assembled matrix (N1 , FI ) 

EM 

Elemental matrix 

N1 

Total number of degrees of freedom 

of the system 

SUBROUTINE BONjDKT 

IE 

Row and column numbers on which 


boundary condition (0) is to be 
applied 
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SUBROUTINE BONDRY 

K Number of boundary conditions 

applied on the system 

N Number of degrees of freed an of the 

system 


SUBROUTINE DYNHAT 


DD 

N1 

Nil 


Dynamical matrix (Nil, Nil) 
Matrix size after applying 
boundary condition 
Size of the dynamical matrix 


SUBROUTINE MATINV 


A 


DETEUM 


Matrix to be inversioned and 
finally its inversion 
Value of the determinant 


SUBROUTINE EIGENP 


A 

EVI(I) 

WR(I) 

N 


Dynamical matrix (NXN) for which 
eigenvalue and eigenvector calculated 
Imaginary part of eigenvalue 
Real part of eigenvalue 
Dimension of the dynamical matrix 



SUBROUTINE EIGENP 


NM 


T 


VECI ( <J,I) 
VECR(J,I) 


Bmensicn ox dynamical matrix given 
in the main program NM >. N 
Number of binany digits in the 
mantissa of a single precision 
floating- point number = 27,0 
Imaginary part of eigenvector 
Real part of eigenvector 


c*r>ooonr*oor*ciooooor*oo 000000000000000 
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APPENDIX - in 


$f,V CUT F OF HP 

$pATCH Z INULA CCSS99 QIO 67 7 

$0UTPU1 U 17 = 1 17 * ’ ' ' 

$ F|\*0F I L 
H LI S Y > 

SIFSYS 
$ I 8FTC MAIN 


FINITE ELEMENT ANALYSIS GF DYNAMIC STABILITY rF PIPE 
CCNvAyING fluid 

DETAILS QF THE PROGRAMME ARE GIVEN IN APPENDIX I AND II 
C AN j ILEyOREQ PIPE HI TF FREE END SUPPORTED 8Y A DISPLACEMENT 
AND A ROTATIONAL SPRING 

fmm is the elemental mass matrix 
tCM IS the elemental camping matrix 

EKM IS THE ELEMENTAL STIFFNESS MATRIX 
AW IS THE OVERALL ASSEMBLED MASS MATRIX 
AC IS THE OVERALL A$SEM8lH0 DAMPING MATRIX 
M IS THF OVERALL ASSEMBLED STIFFNESS MATRIX 
Wl:L IS THE NUMBER QF ELEMENTS 
Or- IS THE DYNAMICAL MATRIX 

ASAMBL IS THE SUBROUTINE' FOR ASSEMBLY QF MATRICES 
V IS THE VELOCITY OF THE FLUID 

oout is the gut side diameter cf the fife 

DIN iS THE INSIDE DIAMETER OF THg PIPE 
8 EL IS THg LENGTH OF THE PIPE 

E IS THE MODULUS OF ELASTICITY CF THg MATERIAL CF TNE PIPE 
Ni-L IS THE NUMBER OF ELEMENTS 

AMF IS THE MASS CF THE FLUID PEr UNIT LENGTH OF THE PIPE 

A M T IS THE MASS OF THg PIPE per unit length cf the pipe 

TA is THE TENSSICN OF THE PIPE INITIALLY 

p IS THE FLUID PRESSURE 

AF IS THE AREA OF THE FLUID PASSAGE 

V'CilS THE nOn-GIHCM signal VELOCITY 

DI IS THE MOMENT OF INERTIA 

ST I FAX IS THE STIFFNESS OF THF DISPLACEMENT SPRING 
ST I F TO IS THE STIFFNESS OF THE ROTATIONAL SPRING 
ALPHA»$TIFAX*(8EL**3 )/<E*Dn=NON-DIMENSIGNAL DISPLACEMENT 
SPRING FACTOR 

BET A=AMF/ (AMF+AMTHNGN-D I m£n SIGNAL MASS 

GAMA»$TIFT0*BEL/(E*D I) “NON-DIMENSIONAL ROTATIONAL SPRING FACTOR 
F “> C N OK ~ ' *OM— D IM EN SI CM A L FAC TOC WHICH WHEN DIVIOEc TO EIGENVALUE 
GIVES NON-DIMENSIONAL FREQUENCY 

FDR CANTILEVER CASE PUT K0=2, IE(1I=1» IE(2J=2 

FOR P IND.-D-P INNfcD CASE PUT KB=2, IE(1)*1, IE C 2 >»NB 


o o o o o o o o o 


IE C 1 ) ~ 1 * IE { 2 1 = 2 


C 


MR FIXED- Fix EC CASE PUT KB=4, 

I ( 3 ) B , I E ( 4 ) = N 1 
F ;,| R I NIL MQ FIXED AND OTHER r ND PI MSEC THEN PUT K8=3» 

!■ < 1 ) = 1 » IE (2 )*2, i E {3)=N6 

WHEN LOTH THE ENDS ARE SUPPORTED BY SPRINGS PUT k 8 = 0 
COMNK,3=OI SpLACEMFNT SPRING AT LEFT E N C 
C ')MNK4 = R0TAT IONAL SPRING AT THE LEFj ChC 
C :iMNK6 = 0I SPIACEMENT SPRING AT THE Rf GFT END 
CGMNM = ROTAT iOKAL spring at the right efo 

0 i H ;.SIUH EMM(4,4) ,ECM4,4) ,EKMl(4,4) ,E (4 , 4 ) , EKH 4,4 ) , AM t 25, 25 ) 

DJMI4SIJN AC(25»25)»Al<(25,25) »0D(5Q*5C) 

DIM'- '.SI Oil FVR{50) ,£VI (50) »VECR(50, 5 O .VECI (50,50 ) * INC IC ( 50 ) 
DIMENSION £VRN(5'* ) ,EVlN(5J» 

DIMENSION IE (1C) 

COMMUN/QYN 1 /AM , AC 
C. -*H /DYN2/D0 
COMMON /DYN 3 /FACTOR 
a: M M CN / D YN 4 / N E l f VNGN 
C-'MMV.^/OYNS/AK 


1 5 


35 


DATA D0UT # DIN»BEl»R0T/0. 024,0.02,2.0,7800.0/ 

DATA E,TA,P/C.2Q4E+1 2,0. 3,0.0/ 

.DATA At PH A, BETA, GAMA /l 0.0, 0.2, 5.0/ 

K J = 2 

PRINT 15, DOy T ,C I N, BE l 

FORM Alt /,5X,*DQUT=*» £15.7, ITERS’* ,2 X«*CIN=*,El5.7 f *M£TERS*,2X, 
l*lj£L=*,el5.7,*M£TERS*l 
PRINT 3 5, TA » P 

FORMAT ( /,5X,*TA**,E15.7»*NEWT0N$*,2X,=«P=* t E15.7,*N/Sq.M*) 
PI=4.0*ATaN(1.) 

CALL FLUN( 20G0Q) 

CALL FLGV( 2GCQC) 

ROF = < OETA/( 1.0-BtTA) ) *RQT* ( (QOlT+DIN I * ( COUT-OIN )/ ( 01 ^**2 ) ) 


PRINT 25,RoTtR0F,E 
FORM AT ( /,5X,*R€T=*,El5.7t 


♦KG/CU.M* ,2X ,*P0F=*,E15.7,*kG/CU.M*, 


22X, n»*»EI5*7**N/SQ«H*) 

AF*Pl*<uiN**2)/4.0 

AT-PI*(D0UT+0IN)«(DoUT-DIN)/4.<1 

AMT «P I*ei;;T* { OCLT+D IN ) * (DOUT-DI M ) /4 .0 

,i' ? F -DCF^AF 

PM:*; !45,AMT,AMF 

45 Fi;TMAT(/,SX»*AMT=# f ElS.7f5X, s *AMF* ! *,E15.7) 

0;-P :*( (DqUT** 2) + (0I N**2 ) )* (DQLT+OIN) *t CCUT-CIN J/64.0 
5uT=SQRT(E*0l/AMFl/BEl 

f ACNH\-$QRTCC6*DI)/< (AMF+AMT)*(BEL**4)) J 
ST IFAX 5 *ALPHA*£*Q I / (861**3 ) 

ST I F7G*GAMA*E*0I /BEL 

PRINT 105*01 » FAC NON, SlIFAXfSTIFTO 

t ,i5 FORMAT(/,5X t *DI-4.61S.7,2X,4FACNON-**6i5.7*2X, 

3$ ST j FAX •♦*61 5* 7,2X »*5T IF 70** f El 5*7) 


n o o 


PRINT 175,AlPHA,£ETA t GAMA 

17 ti I-'GRMaTC /» 5X,*ALPhA=* » E » 7 1 5 X » * 8 ST A=« ,E 15. 7 , 5 X , *G AMA«*. £ 15 *? » 

WRI T :. ( 6» 8$ ) 

86 FORMAT! /,2X»125( 1H*j J 

N !: L a >j 

00 ICO NNEL»1,€ 

£N*N: L 
EL =»iJEL /EM 
PR | NT 135, N6L 
1 3 5 FORmAK /* 2X , #NEL=* »I 3 J 
¥MON=0* 0 
00 30 NV« 1 , l 5 
V V , ":*-SOT 

C'~ r -: \>'~( AMF+AMT )*EL/42C .'0 
CGMNC=aMF*V/ 30.0 
C "'MN ,2*2*C«AMF*V 
CfjMNKl»E*0I/(El**3) 

r'V’O I* ( TA- ( P*AF >-(AMF*(V**2)))/(3G.0*EU 

C0MNK3-STIFAX ->' m 

COMNKAaSTiFlC 

*'••• . ; ®amf*( v**2 )*el 

C/<W'j«C0MNK3 

C0MNK7«CQMNkA 

C0MNK8»(AMF*( V**2) )~TA+(P*AF) 

PR INI' 115, VNCN , V 

lib P">RMaT( /,5X,*VN0N=*,ei5.7»2X,*v a *»615.7 ,*METERS/SEC* ) 

IT (N; L.Nl“.2.0R.VN0N.NE,1.0J GO TO 27 
PRINT 125,C0MNM,C0MNC,C0MNK1,CGMNK2 
125 F »RMAT(/,5X,*COMNM=*,ei5.7,5X,*CciMNC = =«,£15.7,5X,*COMNKl=*,El5.7, 
?5X,*C0MNK2«*,El5.7) 

27 CONTINUE 

Nl*2#(NfcL+l) 

FOrM the ELEMENTAL ma<s matrix 

** * 4 4 si*********************** 3 *** 

fcMMU*n»COMNM*15G.O 
fcMM( 1,2 j sCURNK*22,0*El 
EMM( 1*3)=»C0MNM#54*G 
fcMMi l,4»«-C0)'NM*l>0*Sl. 

(:MM(2,l»*EMMa*2) 

S.-MM ( 2*2 )*C0MNM*4„Q*{ EL**2) 

EMM! 2, 3)»-EMMU»4> 

EMM(2,4)«~CCIMNM*3.0>MEL**2> 

EMM (3,1) “EMM (1,31 
FMM( j,21=EMM(2*3I 
t'^M ( 3,3) “EMM (1,11 
r_«M < j,4 )*-£MMC 1 *2) 

'.MM (4,1) -cMM (1*4) 

IMMI 4,2)»2MM(2*4) 
rW-«, 3)^EMM(3,4) 




nfi o o rr 


, MM ( 4*4 ) =EMM (2*2 ) 

h-rm the Overall mass matrix 
call asahbl(nel, .hm.ni’.amj 

* 'f * *** ‘ ,5t r * * * ™ ^ * ** * * * * * * ** * ** * 3®! * * 3 4 X * 4 * 4 * * 4 # 3( * * * * * 3* # * * # * * * * * * * 

FORM THE ELfcMbNTAL DAMPING MATRIX 

^ ,, ****!! % * < '***** j *** : * ! ** ******** «a*** 

EC.M { i» 1 J =COMNC%3Q.O 
L-CM( 1,2)=CGMNC%6*0*EL 

t-.uM ( 1,3 )=ECM( 1 , 1 ) 

>; CM ( i,4)~~ECM(J. t 2) 

ECM(2,X)=-ECM(l,2) 

LCM( 2*2 »«C*0 
tCMC 2t3}=-6CM(2,i) 

ECM (2,4)=-CGMNC*<EL*#2 ) 

ECM(3,1)*-ECM(1,3) 
cCM (3*2) S -ECM ( 2 1 3 J 
fcCM( 3,31 =-6CH ( 1,1) 
fcCM( 2,4>— ECM3,2) 

ECM< 4 , 1 )*-ECM( 1 , 4 ) 

ECM(4,2)=-ECM(2*4) 

ECM ( <i* 3 1=-ECM( 3*4) 

ECM ( 4 , 4 ) * 0 *G 

c f*jRm the overall damping matrjx 
CALL ASAM8L(NEL,ECM,M,AC) 

AC ( l,U»AC(l.U-COMNC2 
ACl N8» NB ) *AC (NB,NB J+CCMMC2 

C vi „••** ■.’.*&.* ■** ^^-t*«** 4 *$***** * *****4 ******** 4 4 **44 ****** 

C FORM THE ELEMENTAL STIFFNESS MATRIX 

EKMl<i»ll*CQNNKi*n.Q 
EKM H If 2 i *COMNK 1 * 6 * 0 # EL 

F.KMlUt3)»-EKMl( 1,1) 

gKMl( 1,4)®EKM1(1 ,2) 

gKM 1 ( 2* 1 ) =EKM 1 ( 1 ,2 I 

EKM1(2* 2»-C0HNKi*4.0*fEL**2l 

FKM l (2, 3 ) =~EKM1 ( 2*11 

£KM 1 ( 2, 4 ) =COMNK 1*2.0* (6L**2 ) 

£KM1(3»1!--EKM 1(1,3) 

EKMIOf 2)*EKM1(2,3) 

EKM 1 ( 3* 3) *EKMl (1,1) 
fcKM 1 ( 3» A ) *£KM1 (2 ,3 ) 

EKHi(4* 1 ) ®EKM1 (1 ,4) 

EkM 1 ( A« 2 ) “6KM1 (2,4) 

gKM 1(4, 3 )*EKMj(3,4 ) 

£«M 1 ( 4*4,)»£KM1 (2 » 21 

6KM 2 (1*1) *‘C0MNK2*36« 0 
EKM2(1*2)*CQMNK2#3,0*EL 
FKM2(1, 3)»-EK«2(t»U 


non 


C 


c 

c 

c 


c 


, <M2{ 3,4)=?=KM2<lt2> 

iiKM2i 2 t 1 1 =EKM2 ( 1*2) 

fM <M 2 <. 2* 2 ) !EL**2) 

t:KM2!2,3)=-EKM2<2,l) 

t: kM 2 ( 2 » 4 ) =~C GMNK2* (E L*=*2 ) 

j, KM2 ( .3, 1 ) =EKM2 { I » 3 ) 

l:KM,2{ 3f2)=EKM2(2»3) 

*':KM2|3» 3)=cKM2(l*l) 

*.’KM 213,4) =6kM2 (2*3) 

. : <M ;? ( j* * 1 ) = EKM2 { 1, * 4 1 
KM2 ( 4* 2 ) =EKM2 (2 ,4 ) 
i:KM.V (4*3) =fcKM2 (3*43 
? KM2(4,4J=FKM2(2,2) 

0*.' X; 1=1*4 
0 1 1, J = 1 * 4 

EKM( l,J J.EKMHI, J)+EK*2(I t J) 

CONTINUE 

F.3rM THE OVERALL STIFFNESS MATRIX 
CALL ASV'RL (NELtEKM* N1 *AK) 


ak.( i,u=AKa,n+:cK'jK2 
Ak ( 2*2) *AK ( 2*2 )+C 0MNK4 
AK( Ni3»NB J =AK (N8» N6 1+CCMNK6 
AK(N1,NU»AK(N1,NLHCCMNK7 
AK( l»2}=AK(lf2 )-COMNK £ 
AMNB,N1)«Ak<NB»nX>+CCMNK8 


1F(NKL.NE.2.CR.VN0N*NE.1.01 GO TO 37 

WR ITfc! 6* 16 ) 

16 FORMAT! /,12Xt*TH£ OVERALL MASS MATRIX*) 

PR In I 4St ( (AMU* J) t4»l*Nli,l»l,M) 

■»-> FORMAT ( 2X.6E15.7) 

WK ITi. (6*261 

6 format ( /, i2x»*the overall damping matrix*! 

PRINT 55,(<AC(I,J),J = 1,N1),I=1,M) 

5‘j FORMAT ( 2X.6EI5.7! 

WiUTf (6» 36) 

36 FORMAT! /,12X**THE OVERALL STIFFNESS MATRIX*) 

?i f n t 65, (<AK(I, Jl ,j = l,NH»I*ltNl) 

65 FORMAT! 2X,6E15*7) 

.V/ CONTINUE 

APPLYING THE BCuNOARy CONDITIONS 

IS DEFINED EARLIER FOR DIFFERENT BCUACARY CONDITIONS 

N i® 2*(NEL+1 ) 

"t = ht 1 J 

f'sK5=M-KR 

I i' ( K • E Q * 0 ) GO TO 90 


i m=i 
I l 2i=2 
i) • 4 I = I * N l 
4U J “ltNl 
Of .1 ( 1 »J> -AM( f » J) 

a; cunt inul 

call kOndkyc I£»kb*ni > 
n< r : i*i,nkb 

0" 70 J*1 » N«b 

amc i,j)=nou,j) 

70 CONTINUE 

DO 50 1 3 1 »N1 

DO 50 J=1,N1 
DOC I »J)=AC( I » J ) 

50 CONTINUE 

CALL 8QN0RYC IE»K8»N1) 
d ■ an i = i f nkb 
DU a 0 J = 1 » NKB 
ACC 1 1 J I =DO ( I f J) 

' CONTINUE 

DO 60 I ■ I»N 1 
O'! 6 L J 3 1 1 N 1 
DOC I » J ) *AK( I » J ) 

6"> CONTINUE 

CALL nC'l j \Y ( I E ,Kli»Nl ) 

HU 90 1*1, NKB 
DCl 9v J«l,NKB 
AK ( I* J )*D0< I * JI 
■r.,i CONTINUE 

I PC N EL, ME * 2*CFU W: -UN* NE«l* 0 ) GC TO 47 

WRIT; (6,461 „ , 

46 FORMAT* / * 1 2 X * *MASS MATRIX AFTER BOUNDARY CONDITIONS*) 
PRINT 75, ( (ARC I, J),J*1»NKb) t!*I,NK8} 

76 FORM AT < 2X, 2E 15.7 J 

56 fnRMAT?/a2X,*DAMP I NG MATRIX AFTER BOUNDARY CONDITION*! 

PRINT 85,((AC(It J)»4*l »NKBI »I*1 *NKb) 

85 FORMAT* 2X«2t 1.5#7) 

yo { 16*661 

66 f •?. ; ’T{/ f i2X,*$TiFFM*<S MATR iX AFTER EOLNDARY CONDITION*) 
P I '.j 7 95* ( CAKCItOT ,J=1 ,NK8) ,1*1 ,NKBJ 

9b FORHATC 2X » 2E 1 3* 7 ) 

47 CONTINUE 

C FORM THE DYNAMICAL MATRIX 

C ‘ 

CALL DYNMAT(NkB»n1H 

' ‘'if-av * 

•; i. 2 ? » 0 

CALL 6IGENP(Nll,NM,D0iT,6VR f eVl ,YECR, VECI, INCICI 


o n n 


on no i=i f . NU 

0 .-NO*‘l = < (tVRin>»*2) + (:VH 1)4*?J I 

*'VR( I » = cVR( I J/OENOM 

L-vn n=~i~vi( inotiNCM 
;,‘VR'>:m=i r vRcn /facnon 

. V I N ( I J =!; V I ( U/F4CNUN 

1- • CiMTiflU;- 

W 1 T (6,76) 

V ¥' 'KM AT ( /, i 7T f*t:VH*,17> ,*r.VI* ,12X,*II* f lGX,*EVRN*,17X ,*EVIN4 ) 

PRI^T 155, (EVr (I ) tfcVI (I) ,i:VRMl),EVIN(I ) ,1 = 1 » 10 ) 
i5!5 MRMATt llX,£l6.8,3X»EU*B*4X»*II*»4X,E16.tt,5X f E16*8) 

PRINK 165, (ImOIC (I), 1=1, Nil) 

16 5 FORMAT! 5X, 15 1 5 ) 

^ITfc(6,96) 

'■■)h FORMAT! / ,2X,125(1H”) ) 

VNON */*.\T + l.C 
3»’ Ct'iNT INUE 
NEL-N1H. + 1 
WR l Til ( 6 * 106 ) 

106 FORM AT ( /*2X»125(1H#) ) 

100 CONTINUE 
STOP 

END 

$1 BPTC ASAM6L 

Subroutine asambl(nel,em,ni,a«) 

SUBK'uTiNE for THE ASSEMBLY cf matrices 

44444*44444*444*444444**44*44*4***444 44444444444444444444444444444 

DIMENSION EM(4»4}, AM (25,25) 

N"NEL 

M1*2*(N+1J 
00 10 I«1,N1 
DO IQ J*1,N1 
AM! I,J)»0.0 

10 CONTINUE 

DO 2 U INEL«1,N 
N2»(**TN5LI-1 
NT»<2*InEL)+2 
N4»(2*lNEt)-2 
DU 20 I-N2.N3 
00 21 J»N2,N3 
N!»* I-N4 
N6*J-N4 

AMUtJ)*AM(I f J)+EM{N5 ,n6> 

20 CONTINUE 
RETURN 
6nO 

$1 8FTC OYNMAT 

SUOROUTInE OYNHA UNI ,M1) 

r * w jjt #4 *4* ^#4 4 444 44^ 4 44 44 44**#4*f4 44 **♦♦♦* 4 444 


c 

c 

c 

c 

c 


c 

c 


?0 


16 

t~. 5 
3? 


C 

€ 


SUBR.jUT LNE FOR THE ASSEMBLY GF DYNANICAl MATRIX 

* V? * * * 3{s * V * # * * * * * * ♦ t * * * * sft * * * $ * * * * * # * * * :( * * # * * * * * * * * 4 < * * * * * * 4 * * * * * # $ * 

K 18 IS THE REDLCLQ SIZE OF THE MATRIX AFTER gOUNcARY CONCITION 
SCAL il*)G THE STIFFNESS MATRIX AK 

HR i S ‘f HE BIGGEST ELEMENT CF THE STIFFNESS MATRIX AK 
f) IMLNSl ON AM(25,25)»AC(25,25) ,AK(25,2SJ ,00(50 t S0), AMCI25) 
C'JMMCN/QYM 1/AM *AC 
U.V./QYN2/0D 
COMMON/DYN 3/ FACTOR 

0 jM/DYN4/NEL» VNON 
C I N/DYN5/AK 
ftr=abs(ak(i,ih 

O') t<i I =l,Ni 
OU 6 j J-l, Ni 
lF(ABS(AKtI, JJ). 

CONTINUE 
no TO I -1 ,N1 
l)U 70 j * 1 , N i 
AM It-J) =AK ( I , J ) /FTR 
CONTINUE 

COLUNG THE MATRIX INVERSE 
CALL MATINV(fU t DL-TERM) 

1 f ( N r.L. NE • 2 . OR . V NON. NE * l . (J ) 

U}'.; lTfc(6,l6) 

FURmAT(/,12X,*INVERSE of the scaled matfix AK* ) 

PRIM * 25,((AK(!,J),J = I t Ni) ,I = l,Nl) 

FORM AT ( 2X, 2E 15*7 ) 

CONTINUE 

multiply ak Inverse and m and store 


GT.FTP) FTR=ABS(AK(I * J) J 


GO TO 3 7 


jN AM 


OU 2;: 
D'O 30 

SUM*0 
00 40 


J*I,N1 
I ■» l , N 1 

» 0 

K=1,N1 


30 


50 

20 


26 

35 

47 


SUMo SUM+AK (I»K)*AM(K,J) 

CONTINUE 
AHC( I )»SUM/FTR 
CONTINUE 
OU *jO L*1»N1 
AM ( L ,4 ) “AMc ( L ) 

CONTINUE 

continue _ tn 

IF(NEL.NE.2.CR.VN0N.NE.1.0) GO TO 47 
WRITE (6*26) 

p1RMAT</tl2X,*PR0DLCT OF AK INV£R S E AMD AM*) 
PRINT 35, ((AMU, 4) ,J»I»N1) ,1=1, NO 
FORMAT! 2X» 2E1 5.7 j 
CJNTiNUK 

MULTIPLY inverse OF AK AND AC AND STORE IN AC 



u o o 


o'* no j~i»ni 
iji.i no m*Ni 
SOM=0.0 
00 no K=1 , N 1 
SUM = 5Um+AK(I,K)*AC(K, 4> 
l'»n CONTINUE 

AMC l IHSUM/FTR 
120 CONTINUE 

DO HO L«i # NI 
AC ( L , J ) S AMC ( L ) 

14 . CONTINUC- 
11 Oi CO NT t NUE 

IF ( Nir.-L* NE. 2* OR* VfiQ?'!* NE • 1.0 ) 60 Tq 57 
war,. (6,361 

*6 FORMAT! /» 12X ,*PROOUCT OF AK INVERSE Af>D AC* ) 
PRINT 4 5, ( <AC(I,J)»4 = lt*ll),I=l,Nl) 

45 FORMAT! 2X, 2E 15.7} 
hi CONTINUE 

FORM THE DYNAMICAL MATRIX 


C 


BO 


VO 


100 

46 

35 

67 


NlH2*Nl 
DO 6C |*l, Nil 
00 8 C J * 1 ,NI 1 
00 ( I , J } *U *0 
CONTINUE 
Of) SC I»l,Nl 
DO <3U J*1»N1 
Ni'| 1 »N1* I 
NN2°NI+ J 

OOCNNl* J)*-AM( I, J! 

D0( I»NN1H1*C 
DU(NNI,NN2)*-AC( I, J} 

CONTINUE 

FACTOR IS THE FACTqR TAKEN COMMON FROM 

FACTOR* 1.0 


DO ICO I * I , N 1 1 
00 100 J»1,NU 
UD( 1 , J ) *DD ( I , J J/ FACT CP 

CONTINUE 

IF(NfcUN£*2«0R*VNQN»NE*K0J 


GO To 67 


WR I Tf (6,461 

cpiR '-1 AT < /,I2X,*0YNAMICAL MATRIX*) 
PRINT 5 5, <(00(1,4) ,4* J ,NU ) , HI ,Nll ) 
FORMAT ( 2X* 4E 15*7 ) 


CONTINUE 


H TU.V‘) 


THE MATRIX 


M;jO 

IlBpTC MAT INV j _ 

SJ'Hr.UTINE MATINV(N#OE igRM) 






u u u u 


C 

C 

r 

U 

c 

c 

c 


c 


c 


c 


c 


r 


«♦«****« ^ 44 *«:M* 44 ***>|i 44 *# 4 * 444 * 4 *** 

MjOn' ;UT IN'-: FOR THU MATRIX I NVERTION 

'* v. # * ❖ # * # $ # a{c $ t # * ■# * # * 4 ^ ^ 4 s * 4 sjt $ ■# 4 afc * « # * 4 * 4 * ^ => # 4 * $ * * $ * * 4 * j* 4 4 4 4 4 4 4 4 $ 4 4 

kfUTicie FOR MATRIX INVASION a^O SCLUTICN OF SIMULTANEOUS EQN5. 

4 IS THE COEFFICIENT MATRIX CF SIZE U 

0 IS THE VECTOR «jF TH£ SIZE N 

I ; IS NUMBER OF c'^STSM VECTORS 

i 1 M*S SL- T TO ZERO ,QMY INVERSE IS COMPUTED 

O', THEM is the value OF OeTERMENET returned by the routine 

O iMf.ivSIOJl At 25,23) »i PI VUT (40 ) , I NDEX t 4C ,2) 

C0MM0M/DYN5/A 

EQUIVALENCE < IRQW , JRO V ) , ( ICCIUN , JCCLUN ) , (AMAX,T,SWAP ) 
INITIALIZATION 

1 , 0 “ T !f KM= 1 » 0 

15 00 2u J=1,N 

l n I Vl’T ( J ) =0 

SH ARCH FUR PIVOT ELEMENT 
3 D'J 5 30 1 = 1, N 
40 AMAX =CU 0 
4*j 00 U5 Jsl,hl 

5.' 1 1 C IP I VOT < J 1 1 60,105,60 

bo DU ICO K= 1 , N 
7 ,' IK I PI VOT ( K 1-1 ) SO, ICC ,740 
SO II ( AMAX-ABSC A(J,K) n 8 5,100,100 

1 UW - J 

VO I COLUMN 
9 i At:AX=ABS(A( J,KH 

iOO CONTINUE 
105 CONTINUE 

IlC IP I VUT ( ICOLUM) »I PI VOT (ICOLUM J + l 

INTERCHANGE THE ROWS TO PIT PIVOT ELEMENT CN DIAGONAL 
130 I Ft IRQW-ICOLUM) 140,2*0, HO 
140 OETFRM— OETERM 
150 00 2C€ L » 1 * N 
160 SWAP»A( IRON, LI 
170 AURoW»L)»A<ICCLt w ,L 1 
20',' AC ICULUM»LI*SwAP 

260 INOEXCI,U*lROW 

2 ?c rorxu, 2 )»icolu« 

DIVIDE PIVOT Row BY PIVOT SL'-MqNT 
3iC P ? VO T*A ( ICOLUM * ICOLuN J 

329 OFT: CM-0uTFRM*P1V0T 

330 AC ICULUM, IC0LUM1 *1«0 
AO. ui 350 L*l,N 

,4., AC i CULUH*L)*AUCOLIM,U/PIVCt 
RvC’JCC NON PIVOT ROWS 
ibr ,-j 1 550 Ll*l,N 
iv ,1 IF(Ll-ICOLUM) 400,550,400 
•’> k! T =A ( L l , 1C0LUM ) 


no o o o o n o n o 


4 }' f, A(U,lCOLUM|r(}.Q 

43*' 0, J 431,* 1=1, N 

4 5!.' A(U,U=A(Ll t U-A(ICCtUN,U*T 

55'' C NT li'iUl: 

C r NT 1; ".CHANG COLUMNS 

6tV1 0:i 710 1 = 1, N 
610 L=N+1-| 

620 IF{ «Di:X(L,l)--lNUEX( 1,21) 630,710,630 
6 3:. ji-:ow = iNOs. : xa,i) 

6'*' !') 4C0LliM*IN0EX (L ,2 J 

650 mo 705 K=L,N 

660 5 sAP=A(K, JR04O 

t'7C A ( K , JRO W ) = A ( K» JCCLUM > 

7 1 • C A (K , JCOLUM ) = SWAP 
705 CJ\T1NU6 ' 

/ 10 CONTINUE 

DO 11 K = 1 , N 

l F ( IPiVOT(K).NE.n GO TO 12 

Li CONTINUE 
K, TURN 

12 PRINT 991 

991 H!RMAT< /30X*MATR1X IS S* ?i| CCL *?*/) 

7*v0 '< • ' T y 
t-iO 

$IBT T C BON DRY 

SUBROUTINE 80N0RY( 16 ,K,N) 

** 4 ^ 4 $ *$$**#* $#**$* 

> 3“ I’.IT;..; FOR THE BOUNDARY CONDITION 

jjt 4 ** $* 4 c****i|c *jM##*$#** ***** > 4< -$<*#*&* **4 taM************************** 

Ism IS THE ROW OR COLUMN NUMBER FOR WHCH THE 
BOUNDARY CONDITION IS APPLIED AND IS 2£PC 
K IS THE TOTAL NO. OF BOUNDARY CONDITION APPLIEC 
N IS THE SIZE OF THE ORIGINAL MATRIX EEFCFE APPLYING THE 
BOUNDARY CONDI 1I0NS 
00 IS THE MATRIX 

D I M N S I ON 16(41,00(50*50) 

r, mm ;n/dyn2/oo 
kr«c 

00 1(1 1 = 1 , N 
UU 2C 1*1, K 

IFC I.fcQ.IE(L) ) GO TO 10 

20 Continue 

KR«KR+1 
kC w G 

00 4C 4=1 , H 
00 3C l-1,K 

IF( J.EQ.IE(L ) ) GO TO 40 
30 CONTINUE 
KC*KC+1 




i;d(kr,kc)=cd<i,j) 

*:• C'.IMT i'NUii 
i C 1 NT HU if 


»v TLH 

r o 

& 1 U, i.j L T 

tlulin :: I G . K‘P 
J13LDP Gc A L )~ 

$ i uL OR HESWR 
ii^tCR '^alve 
S £ DL. OR COMPVE 
tc'ir^y 
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